2日目
SPSSでは,できない解析がある
機能のわりに高価(卒業後は使えない)
操作マニュアルは多いが,教育上は悪い影響も(統計学を理解せずに解析するなど)
データの整形機能が貧弱
解析過程が残らない(解析が再現出来ない可能性がある)
出力が非常に見にくい(好みの問題ですが・・・)
SPSSでデータ整形は難しいので,Excelで前処理をすることも多いのでは?
しかし,以下の理由から,Excelもすすめない。
操作の過程が残らない(Excelで行った操作を別のファイルに記録しないと再現不可能になる)
生データ自体が改変される
本格的な統計解析・可視化は難しい(HADという選択もあるが・・・)
できない統計解析はほとんどない
解析過程がしっかり残る
データの整形,可視化,解析,解析レポート作成(マニアは論文やプレゼンも)可能
こんなにすごいのに,無料!
学生が統計を学習する上で便利(自宅でもできる&修了後も使える)
The R Project for Statistical Computingにアクセスし,「download R」をクリックする。
Rtoolsのインストール
Windowsユーザーの注意事項
ユーザー名を日本語名にしていると問題が発生することがある。英語名のユーザーを追加するか,こちらを参考にして,設定を変更する。
Xcodeのインストール
コマンドラインツールのインストール
xcode-select --install
Xquartzのインストール
Rのインストール
Rは標準関数で色々な統計解析ができる。
Positチートシートの「Base R」(日本語版)を参照したり,「R 使いたい解析名」などでグーグル検索して調べる。
Rの関数名はわかっているけど,どう使ったらわからない場合,「?関数名」をコンソールにタイプするとヘルプがでてくる。
戻り値(結果) = 関数(引数(入力))
例) 平均値の関数mean()の場合
結果 = mean(平均したいデータ)
Rでは,標準関数以外にも様々な解析パッケージが作成され,多くはCRANに登録されている。
CRANにあるパッケージは以下のようにコンソールにタイプするとインストールできます。
install.packages("パッケージ名", dependencies = TRUE)
library(パッケージ名)
今回の講義で使用するパッケージをインストールする。
dmetarはgithub経由でインストールする。
パッケージを探す:「R 解析名」でグーグル検索
関数の使い方を探す:(1)「?関数名」をRStudioのコンソールにタイプするとヘルプがでてくるので読む,(2)CRANでパッケージを検索してマニュアルを読む(分かりやすいビネットが公開されていることもある),(3)「関数名 解析名」などでグーグル検索(英語のほうが情報が多そうだがRは日本語でも重要な情報が見つかることも多い)
Consoleにコマンドを打ち込むやり方だと手元に解析過程が残らない,明示的にコマンドとその説明文章を残す習慣をつける必要がある。
Quartoを使うと,便利なMarkdown記法を使って,文章と解析コードを残せる。
Quartoは複数の出力(HTML, PDF, Word)ができ,R以外の言語も対応し,RStudio以外でも利用できる(Quartoは元々あったRmarkdownを統合したものという位置づけ)。
SourceモードとVisualモードがあるが,VisualモードはWordみたいな感じで作業ができる。
①設定情報。②本文(見出しや引用や画像の挿入)。③解析部分のチャンク。④RenderをクリックするとHTML出力される。
RとRStudioとRとパッケージをインストールして,Quartoドキュメントを作成して(書き込みはしなくてもOKです),Renderしてみよう!
時間が空いたなら,以下の資料を読んで再現可能な環境構築について学ぼう。
統計学やプログラミングは実際に手を動かして学ぶのが効率的
Rのパッケージには統計の学習用のデータも用意されていますが,キレイ整いすぎている。そこで,実際のオープンデータを使うのをオススメする。オープンデータを使った解析は,2次分析研究にもなり得る。
2次分析(Secondary Data Analysis)とは,既存データ(preexisting data)に,既存の研究とは違う目的(視点)をもって分析することである。
大規模調査のデータ(政府統計,パネル調査)
個別の研究データ(各研究室のデータ)
ビッグデータ(ウェブスクレイピング,医療報酬の明細書など)
既存データの中にはオープンデータも増えている。
オープンデータとは,「自由に使えるデータ」(庄司, 2014)
→どのようなライセンスの下でデータが公開されているのかが重要になる
オープンデータのリストとデータ検索サイト
データを掲載するジャーナルも増えてきている。
国里も参加した新型コロナ禍中の国際共同プロジェクトのデータはScientific Dataに掲載(世界で73,223名が参加した研究)
対象者:一般参加者295名
介入:3つのポジティブ心理学的介入と統制条件に無作為割付。(1)“Using Signature Strengths(ウェブサイトで性格の強みを評価し,上位の強みのうちの1つを次の週に新しい別のやり方で実施する)”,(2)“Three Good Things(その日に起こった良いことを3つあげて,なぜそれが起こったのか原因説明とともに記録する)”, (3)“Gratitude Visit(親切してくれたけど感謝を示してない人に手紙を書く介入)”,(4)“Recording early memories(統制条件,幼いことを記憶を書く)”
アウトカム:Authentic Happiness Inventory (AHI) の24項目,Center for Epidemiological Studies Depression(CES-D) scaleの20項目
測定時期:介入前(0), 介入後(1), 1週間後(2), 1ヶ月後(3),3ヶ月後(4),6ヶ月後(5)
Woodworth et al.(2018)のデータにアクセスして,ファイルをダウンロードして開きます。
RStudioにはプロジェクトという機能があり,フォルダやファイルやパッケージの管理に便利です。研究プロジェクト単位でプロジェクトを作って,関連するファイルをいれると解析に便利です。
FileからNew Projectをクリック
設定部分より下は削除して,日本語入力してみましょう。書式設定とかできるので,適当にいじってみましょう。また,見出しはとても重要です。文章の構造化を意識しましょう。
Rチャンクの追加は,Insert -> Code Cell -> Rでできます。
データの読み込み,データの前処理,可視化などはtidyverseパッケージ使うと便利です。
dataフォルダに保存したWoodworth et al.(2018)のデータを読み込みます。データはcsv形式だったので,read_csv()で読み込みます。ファイルはparticipant-info.csvとahi-cesd.csvの2種類があります。dataフォルダにいれているので,ファイル名の前に”data/“をいれています。こうすると,dataフォルダ内の特定のファイルを指定できます。
Rではパイプ演算子という便利なものを使います。以前はmagrittrパッケージを用いた%>%だったのですが,R4.1からは標準で|> が使えます。
|> が使えるように設定します。Tools -> Global options -> Codeで,Use native pipe operator…にチェックをいれてApplyする。
Windowsの場合Ctrl + Shift + m Macの場合Command + Shift + mのショートカットでパイプ演算子が入力できます。
Rではデータフレームという構造でデータを扱う。データフレームは,a行b列からなる表形式のデータであり,列には変数が含まれる。
read_csv()で読み込んだデータはデータフレームになっており,以下のオブジェクト名を打つとデータフレームが表示される。
| データ型 | 説明 | 省略形 |
|---|---|---|
| logical | 論理型(TRUE, FALSE) | lgl |
| integer | 整数型(1,2) | int |
| double | 倍精度小数点型 (1.0,2.1) | dbl |
| character | 文字型 | chr |
| factor | 因子型 | fct |
| ordered | 順序型 | ord |
| Date | 日付型 | date |
前処理にはtidyverseパッケージに含まれているdplyrパッケージを使う。色々と機能はあるが,ミニマムには以下の8つを覚える。
filter() 行の絞り込み
arrange() 行の並び替え
select() 指定した列の抽出
mutate() 列の追加
group_by() グループ化
summarize() 要約
pivot_longer() 横データから縦データへ変換(時系列解析などは縦データが便利)
pivot_wider() 縦データから横データへ変換
filter() で行の絞り込みができます。性別を女性(1,男性は2)だけに絞り込む場合は以下のようにできます。パイプ演算子は,左にあるオブジェクトを右(改行していることが多い)の関数にいれる機能を持っています。パイプ演算子を使うと処理を重ねていけるので見やすいし便利です。今回は,participantデータをfilter()にいれて,sexが1のものだけを抽出しています。
パイプ演算子の挿入はWindowsの場合Ctrl + Shift + m Macの場合Command + Shift + mです。
# A tibble: 251 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2 1 1 59 1 1
2 3 4 1 51 4 3
3 4 3 1 50 5 2
4 6 1 1 31 5 1
5 7 3 1 44 5 2
6 8 2 1 57 4 2
7 9 1 1 36 4 3
8 10 2 1 45 4 3
9 11 2 1 56 5 1
10 12 2 1 46 4 3
# ℹ 241 more rows
summarize() を使うと,変数名=関数(変数)で要約ができます。最小40歳,最大85歳ですね。# A tibble: 1 × 2
min_age max_age
<dbl> <dbl>
1 40 83
以下のように40歳以下かつ女性で絞り込むこともできます。
# A tibble: 17 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 4 2 35 5 3
2 15 4 2 27 2 2
3 44 2 2 36 4 3
4 130 4 2 21 1 1
5 140 3 2 31 3 1
6 193 2 2 25 5 2
7 204 4 2 32 5 1
8 210 1 2 28 4 2
9 234 3 2 33 5 3
10 235 2 2 18 2 1
11 237 4 2 33 3 3
12 243 4 2 26 4 2
13 251 4 2 31 5 2
14 263 1 2 32 5 2
15 269 3 2 23 4 1
16 284 3 2 19 4 1
17 288 3 2 26 5 1
arrange() で行の並び替えができます。指定した変数で昇順で並び替えます。# A tibble: 295 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 235 2 2 18 2 1
2 267 4 1 19 2 1
3 284 3 2 19 4 1
4 20 3 1 21 4 1
5 130 4 2 21 1 1
6 230 2 1 21 2 1
7 23 1 1 22 4 1
8 122 1 1 23 5 2
9 138 1 1 23 2 1
10 247 1 1 23 4 1
# ℹ 285 more rows
# A tibble: 295 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 51 2 1 83 2 2
2 244 2 1 77 3 2
3 215 4 1 75 3 2
4 83 4 1 71 2 1
5 114 1 1 71 4 1
6 155 3 2 71 4 1
7 93 3 1 68 1 2
8 141 2 2 67 5 2
9 281 4 2 67 4 2
10 67 1 1 65 4 2
# ℹ 285 more rows
select()で指定した列の抽出ができます。例えば,id,age,sexだけ抽出したい場合は,以下のように選択できます。mutate() で列の追加ができます。例えば,dataにはCES-Dの得点が入っていますが,その中でポジティブ感情の項目(4,8,12,16)だけを合計して得点化します。確認のためにsummarizeもしています。group_by()でグループ化します。dataのCES-Dの合計得点(cesdTotal)の介入前の得点(occasionが0)に対して,interventionでグループ分けして,summarizeします。pivot_longer()で横データから縦データへ変換,pivot_wider()で縦データから横データへ変換します。dataはすでにlong形式になっています(occasionごとに横にデータが並んでいるのではなく,1つの参加者でoccationごとに複数行が並んでいます)。これを横データにしてみます。select(id,occasion,intervention,cesdTotal)で少し変数をしぼって,pivot_wider(names_from = occasion, values_from = cesdTotal) で横データにします。最後に,-> data_wideとするとdata_wideというオブジェクト名として保存されます。occasionの1以降は回答してない人も多いので,NULLになっていますね。
idの8がoccasionの2で二回回答しているので,二回目の回答を除外します(指定がしにくくてid とcesdTotalで指定しています)。
pivot_longer()で縦データを横データに変換します。pivot_longer(cols = c("移動する列"), names_to = "分ける変数名", values_to = "値の名前") で横データになります。ただ,元に戻るのではなく,横→縦変換時に欠測になっているので,横データに戻した時に欠測が含まれています。data_wide |>
pivot_longer(cols = c(-id,-intervention), names_to = "occasion", values_to = "cesdTotal")# A tibble: 1,770 × 4
id intervention occasion cesdTotal
<dbl> <dbl> <chr> <list>
1 1 4 0 <dbl [1]>
2 1 4 1 <dbl [1]>
3 1 4 2 <NULL>
4 1 4 3 <NULL>
5 1 4 4 <NULL>
6 1 4 5 <NULL>
7 2 1 0 <dbl [1]>
8 2 1 1 <dbl [1]>
9 2 1 2 <dbl [1]>
10 2 1 3 <dbl [1]>
# ℹ 1,760 more rows
Woodworth et al.(2018)のデータについて,8つの関数を駆使して,様々にデータを変形してみましょう!
Rows: 295
Columns: 6
$ id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ intervention <dbl> 4, 1, 4, 3, 2, 1, 3, 2, 1, 2, 2, 2, 4, 4, 4, 4, 3, 2, 1, …
$ sex <dbl> 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, …
$ age <dbl> 35, 59, 51, 50, 58, 31, 44, 57, 36, 45, 56, 46, 34, 41, 2…
$ educ <dbl> 5, 1, 4, 5, 5, 5, 5, 4, 4, 4, 5, 4, 5, 1, 2, 1, 4, 5, 3, …
$ income <dbl> 3, 1, 3, 2, 2, 1, 2, 2, 3, 3, 1, 3, 3, 2, 2, 1, 2, 2, 1, …
gtsummary パッケージのtbl_summary()を使うと簡単に変数の要約ができます。| Characteristic | N = 2951 |
|---|---|
| intervention | |
| 1 | 72 (24%) |
| 2 | 76 (26%) |
| 3 | 74 (25%) |
| 4 | 73 (25%) |
| sex | |
| 1 | 251 (85%) |
| 2 | 44 (15%) |
| age | 44 (34, 52) |
| educ | |
| 1 | 14 (4.7%) |
| 2 | 21 (7.1%) |
| 3 | 39 (13%) |
| 4 | 104 (35%) |
| 5 | 117 (40%) |
| income | |
| 1 | 73 (25%) |
| 2 | 136 (46%) |
| 3 | 86 (29%) |
| 1 n (%); Median (IQR) | |
byでグループ分けをして表を作ることもできる。| Characteristic | 1, N = 721 | 2, N = 761 | 3, N = 741 | 4, N = 731 |
|---|---|---|---|---|
| sex | ||||
| 1 | 62 (86%) | 66 (87%) | 62 (84%) | 61 (84%) |
| 2 | 10 (14%) | 10 (13%) | 12 (16%) | 12 (16%) |
| age | 45 (36, 54) | 46 (36, 53) | 44 (34, 51) | 40 (32, 51) |
| educ | ||||
| 1 | 3 (4.2%) | 3 (3.9%) | 2 (2.7%) | 6 (8.2%) |
| 2 | 3 (4.2%) | 8 (11%) | 3 (4.1%) | 7 (9.6%) |
| 3 | 10 (14%) | 7 (9.2%) | 7 (9.5%) | 15 (21%) |
| 4 | 26 (36%) | 32 (42%) | 32 (43%) | 14 (19%) |
| 5 | 30 (42%) | 26 (34%) | 30 (41%) | 31 (42%) |
| income | ||||
| 1 | 21 (29%) | 20 (26%) | 17 (23%) | 15 (21%) |
| 2 | 31 (43%) | 34 (45%) | 36 (49%) | 35 (48%) |
| 3 | 20 (28%) | 22 (29%) | 21 (28%) | 23 (32%) |
| 1 n (%); Median (IQR) | ||||
gtsummary()の調整法としては,こちらが参考になります。
より柔軟に表を作成するパッケージとしてはgtパッケージがありますが,説明を省略します。
2群の差は,標準化平均値差(Hedgesのg)を用いる
\(Hedgesのg = \frac{介入群の平均-統制群の平均}{2群をプールした標準偏差}\)
\(2群をプールした標準偏差=\sqrt{\frac{(n_{介入}-1)\hat{\sigma}^{2}_{介入}+(n_{統制}-1)\hat{\sigma}^{2}_{統制}}{(n_{介入}+n_{統制}-2)}}\)
\(\hat{\sigma}^{2}\)は不偏分散。他にCohenのdやGlassの⊿もあるが,バイアスを修正した効果量のHedgesのgがよく使用される
\(標準化平均値差の分散=\frac{n_{介入}+n_{統制}}{n_{介入}n_{統制}}+\frac{g^2}{2(n_{介入}+n_{統制})}\)
\(標準化平均値差の標準誤差=\sqrt{標準化平均値差の分散}\)
分散と標準偏差は効果量とサンプルサイズで決まる
gtパッケージのgt()できれいにした。library(effsize)
library(gt)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
select(intervention,cesdTotal) |>
group_by(intervention) |>
summarize(mean = mean(cesdTotal),sd = sd(cesdTotal), n = n()) |>
gt()| intervention | mean | sd | n |
|---|---|---|---|
| 3 | 16.08108 | 11.815658 | 74 |
| 4 | 12.84932 | 9.511209 | 73 |
Hedgesのgを計算すると弱い効果量がある。
Cohen's d
d estimate: 0.3010943 (small)
95 percent confidence interval:
lower upper
-0.0267866 0.6289752
2変数間の関連の強さは,ピアソンの積率相関係数を用いる
\(r = \frac{共分散_{x,y}}{標準偏差_{x}標準偏差_{y}}\)
\(相関係数の分散 = \frac{(1-r^2)^2}{n-1}\),\(相関係数のSE = \frac{1-r^2}{\sqrt{n-1}}\)
統計学的特性の望ましさから(サンプルサイズが小さいとrのseにはバイアスが生じる),フィッシャーのZ変換を用いることが多い
\(z = 0.5ln(\frac{1+r}{1-r})\), \(zの分散 = \frac{1}{n-3}\), \(zのSE = \frac{1}{\sqrt{n-3}}\)
zからrへは\(\frac{exp(2z)-1}{exp(2z)+1}\)で変換できる(信頼区間も同様)。
リスク差=介入群のありの比率-統制群のありの比率
\(治療必要数=\frac{1}{リスク差}\)
\(リスク比=\frac{介入群のありの比率}{統制群のありの比率}\)
\(オッズ比=\frac{\frac{介入群のありの比率}{1-介入群のありの比率}}{\frac{統制群のありの比率}{1-統制群のありの比率}}\)
リスク差や治療必要数(Number needed to treat)は解釈しやすいが,計算上オッズ比が便利なので,オッズが使用される。さらに対数オッズ比(ln(オッズ比))を使う(1を起点に解釈→0を起点に解釈,分散も計算がしやすい)。
オッズ比の分散=介入群あり比率+介入群なし比率+統制群あり比率+統制群なし比率
Woodworth et al.(2018)のデータについて,論文のTable1に該当するような表を作成してみよう!他の群間の平均値を比べて,効果量を計算してみよう!
Rには可視化のための関数があるが,ggplot2パッケージがきれいに図が作成できて,コードも読みやすく書ける。
ggplot2は,グラフィックの文法 (grammar of graphics)という思想の元に作成されている。これは,有る決まった構造で作成されたレイヤーを重ねていって図を作成するという考え方です。
+を使います。aes()では,x軸, y軸の変数や色などを指定します。ここでは,1変数のsex2だけ指定します。macの場合は文字化けしているのですが,ちょっと置いておきます。geom_bar()で,棒グラフを重ねます。base_family = "HiraKakuPro-W3" を指定します。theme_classic()を使います。あれ?グラフが浮いている感じして気持ち悪いですね。ggplot(data = participant) +
aes(x=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") scale_y_continuous(expand=c(0,0)) を追加すると浮かなくなります。ggplot(data = participant) +
aes(x=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0))aes()内のfillにsex2を指定します。ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0))ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray"))あ!線がない!geom_bar()内に色を追加します。
ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar(color="black") +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray"))labs()を追加します。ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar(color="black") +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray")) +
labs(x = "性別", y = "人数")theme(legend.position = "none")を追加します。ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar(color="black") +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray")) +
labs(x = "性別", y = "人数") +
theme(legend.position = "none")使えそうな図が作れました。ggplot2は面倒な感じもしなくもないですが,キャンバスにちょっとずつレイヤーを足していくことで図を作る感じをイメージするといいです。その都度,「ggplot 凡例 消す」とかで検索すると簡単に情報が手に入ると思います(chatGPTも便利です)。
量的な変数についてはヒストグラムを書いて分布を確認します。
dataからoccasion=0のCES-D得点のヒストグラムを描きます。filterでoccasion=0に絞り込んで,それをggplotに投げます。パイプ演算子を使っているので,filterされたdataがggplotに入るので,aesだけ指定します。geom_histogram()でヒストグラムが描けます。
data |>
filter(occasion == 0) |>
ggplot(aes(x=cesdTotal)) +
geom_histogram() +
theme_classic() +
scale_y_continuous(expand=c(0,0)) +
labs(x = "CES-D") geom_histogram(aes(y=after_stat(density)),fill='white', color='black') という指定をしています。また,geom_density(fill='black', alpha = 0.5) のalphaは重ねた際の透過度になります。geom_point()で描けます。data |>
filter(occasion == 0) |>
ggplot(aes(x = cesdTotal, y = ahiTotal)) +
geom_point() +
theme_classic() +
labs(x = "CES-D", y = "AHI")data |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
ggplot(aes(x = cesdTotal, y = ahiTotal,
fill = intervention, color = intervention)) +
geom_point() +
theme_classic() +
labs(x = "CES-D", y = "AHI")時間的な変化や回帰直線などをプロットする際には,折れ線が使える。
今回は,occasionが0~5まで5時点あるので,id=1〜10までプロットしてみる。工夫としては,idで10以下のデータにしぼったこと,idは因子型にしたことになる。折れ線は,geom_line()を使う。
data |>
filter(id <= 10) |>
mutate(id = as.factor(id)) |>
ggplot(aes(x = occasion, y = cesdTotal, color = id)) +
geom_line() theme_classic()などの見栄えの調整をした。geom_point()も重ねると点がわかりやすいので追加している。geom_boxplot()でプロットできる(stat_boxplot()はエラーバーを出すため)。data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
ggplot(aes(x = intervention, y = cesdTotal)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot()data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot() +
theme_classic() +
labs(x = "Intervention", y = "CES-D")geom_violin(trim = FALSE)で描けます。data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
geom_violin(trim = FALSE) +
theme_classic() +
labs(x = "Intervention", y = "CES-D")data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
geom_violin(trim = FALSE) +
geom_jitter(position=position_jitter(0.1))+
theme_classic() +
labs(x = "Intervention", y = "CES-D")stat_halfeye()でバイオリンプロットの半分みたいなものを描きます。library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye()stat_halfeye()の信頼区間っぽいやつを削除するために,point_color = NA, .width = 0を指定します。library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.05, binwidth = 1)library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
geom_boxplot(width = .1,outlier.shape = NA) coord_cartesian(clip = 'off', xlim = c(1, 2))で調整しました。library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
geom_boxplot(width = .1,outlier.shape = NA) +
theme_classic() +
labs(x = "Intervention", y = "CES-D") +
coord_cartesian(clip = 'off', xlim = c(1, 2))case_when()を使っています。library(ggdist)
data |>
filter(occasion == 0) |>
mutate(intervention2 = case_when(
intervention == 1 ~ "Signature Strengths",
intervention == 2 ~ "Three Good Things",
intervention == 3 ~ "Gratitude Visit",
intervention == 4 ~ "Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
geom_boxplot(width = .1,outlier.shape = NA) +
theme_classic() +
labs(x = "Intervention", y = "CES-D") +
coord_cartesian(clip = 'off', xlim = c(1, 4))データの分布が正規分布:t検定
サンプルサイズ小さい,データの分布が正規分布から大きく外れる:ウィルコクソン順位和検定
データの分布が正規分布の場合t検定を使う。
Rのt検定では,t.test()を使う。
まずは平均と標準偏差を並べる。
| intervention | mean | sd |
|---|---|---|
| 3 | 16.08108 | 11.815658 |
| 4 | 12.84932 | 9.511209 |
t検定は,t.test(従属変数~独立変数, data = data)で実施できるが,以下の引数がある。
alternative = c(“two.sided”, “less”, “greater”) : 基本は両側検定
mu = 0:もし1標本t検定を行う場合
paired = FALSE:対応のあるt検定を行う場合
var.equal = FALSE:両群の等分散性が確認できないと通常のt検定はできず, Welch の補正 (Welch’s correlation)が必要になる。rではデフォルトでWelch の補正をかけている。
まず,必要なデータを整理して,data_2sample_tに保存します。際どいですが,有意差はないようです。ただ,無視できない差に思います。
data_2sample <- data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
select(intervention, cesdTotal)
t.test(cesdTotal ~ intervention,
data = data_2sample,
alternative = "two.sided",
paired = FALSE,
var.equal = FALSE)
Welch Two Sample t-test
data: cesdTotal by intervention
t = 1.8279, df = 139.41, p-value = 0.0697
alternative hypothesis: true difference in means between group 3 and group 4 is not equal to 0
95 percent confidence interval:
-0.263802 6.727334
sample estimates:
mean in group 3 mean in group 4
16.08108 12.84932
データに正規性に疑いがある場合,ノンパラメトリック検定のウィルコクソンの順位和検定を用いる。ウィルコクソンの順位和検定では,2群の順位から検定を行う。
Rのウィルコクソンの順位和検定では,coinパッケージのwilcox_test()を使う。
Exact Wilcoxon-Mann-Whitney Test
data: cesdTotal by intervention (3, 4)
Z = 1.4016, p-value = 0.1618
alternative hypothesis: true mu is not equal to 0
2群の比率の検定では,フィッシャーの正確検定を行います。
まずはクロス集計を計算します。
data_2ratio <- participant |>
filter(intervention == 3 | intervention == 4)
table(data_2ratio$sex, data_2ratio$intervention)
3 4
1 62 61
2 12 12
fisher.test()で実行できます。
Fisher's Exact Test for Count Data
data: data_2ratio$sex and data_2ratio$intervention
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.3841867 2.6885100
sample estimates:
odds ratio
1.016284
chisq.test()で実行できます。サンプルサイズが十分に大きい,もしくはデータが正規分布:1要因分散分析
サンプルサイズが小さい,データが正規分布ではない:クラスカル・ウォリス検定
もし1要因分散分析やクラスカル・ウォリス検定で有意差があれば,多重比較を行う。
1要因分散分析(一元配置分散分析)は,oneway.test()で実行できる。
dataのcesdTotalを4群で比較する。まずはデータを用意して,平均と標準偏差を確認する。
data_4group <- data |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
select(intervention, cesdTotal)
data_4group |>
group_by(intervention) |>
summarize(mean = mean(cesdTotal),
sd = sd(cesdTotal)) |>
gt()| intervention | mean | sd |
|---|---|---|
| 1 | 15.05556 | 11.677460 |
| 2 | 16.21053 | 9.894532 |
| 3 | 16.08108 | 11.815658 |
| 4 | 12.84932 | 9.511209 |
data |>
filter(occasion == 0) |>
mutate(intervention2 = case_when(
intervention == 1 ~ "Signature Strengths",
intervention == 2 ~ "Three Good Things",
intervention == 3 ~ "Gratitude Visit",
intervention == 4 ~ "Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
geom_boxplot(width = .1,outlier.shape = NA) +
theme_classic() +
labs(x = "Intervention", y = "CES-D") +
coord_cartesian(clip = 'off', xlim = c(1, 4))oneway.test()で決定を行います。
One-way analysis of means (not assuming equal variances)
data: cesdTotal and intervention
F = 1.8107, num df = 3.00, denom df = 160.69, p-value = 0.1473
kruskal_test()で実行できる。
Asymptotic Kruskal-Wallis Test
data: cesdTotal by intervention (1, 2, 3, 4)
chi-squared = 4.6244, df = 3, p-value = 0.2015
pairwise.wilcox.test(data_4group$cesdTotal, data_4group$intervention, p.adjust.method = 'bonferroni')
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: data_4group$cesdTotal and data_4group$intervention
1 2 3
2 1.00 - -
3 1.00 1.00 -
4 1.00 0.19 0.97
P value adjustment method: bonferroni
3群以上の比率の差は,カイ二乗検定を用いる。カイ二乗検定は,chisq.test()を用いる。
dataの介入群ごとの性別の比率をクロス集計で確認する。
Pearson's Chi-squared test
data: participant$intervention and participant$sex
X-squared = 0.47685, df = 3, p-value = 0.9239
pairwise.prop.test()で行います。2つの連続変数の関係は散布図や相関係数で表現できる。それが無相関かどうかの帰無仮説検定するのが無相関検定である。
まずは散布図を確認する。
data |>
filter(occasion == 0) |>
ggplot(aes(x = cesdTotal, y = ahiTotal)) +
geom_point() +
theme_classic() +
labs(x = "CES-D", y = "AHI")cor()で計算できる。相関係数は2つの連続変数間関係についての効果量である。
Pearson's product-moment correlation
data: data_cor$ahiTotal and data_cor$cesdTotal
t = -18.042, df = 293, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.7754145 -0.6664727
sample estimates:
cor
-0.7254563
基本的にはcor()やcor.test()で十分だが,多変数の相関関係を確認したい場合は,corrrパッケージが便利かもしれない。
CES-Dの20項目の相関関係を検討する。correlate()で相関行列を計算して,cesd_corに格納する。
fashion()で見やすくできる。 term cesd01 cesd02 cesd03 cesd04 cesd05 cesd06 cesd07 cesd08 cesd09 cesd10
1 cesd01 .38 .43 -.21 .41 .39 .37 -.17 .27 .29
2 cesd02 .38 .40 -.23 .20 .35 .29 -.18 .31 .18
3 cesd03 .43 .40 -.44 .43 .73 .59 -.42 .52 .44
4 cesd04 -.21 -.23 -.44 -.26 -.40 -.32 .57 -.33 -.23
5 cesd05 .41 .20 .43 -.26 .45 .46 -.22 .27 .29
6 cesd06 .39 .35 .73 -.40 .45 .67 -.45 .59 .49
7 cesd07 .37 .29 .59 -.32 .46 .67 -.31 .49 .42
8 cesd08 -.17 -.18 -.42 .57 -.22 -.45 -.31 -.40 -.32
9 cesd09 .27 .31 .52 -.33 .27 .59 .49 -.40 .45
10 cesd10 .29 .18 .44 -.23 .29 .49 .42 -.32 .45
11 cesd11 .28 .25 .34 -.18 .21 .37 .30 -.19 .21 .30
12 cesd12 -.41 -.35 -.64 .49 -.38 -.65 -.53 .56 -.43 -.44
13 cesd13 .33 .28 .40 -.23 .29 .35 .29 -.21 .26 .33
14 cesd14 .26 .31 .45 -.32 .26 .49 .39 -.34 .52 .38
15 cesd15 .23 .09 .16 -.09 .10 .16 .21 -.07 .28 .17
16 cesd16 -.31 -.33 -.60 .48 -.31 -.56 -.43 .53 -.40 -.35
17 cesd17 .27 .29 .49 -.18 .30 .36 .37 -.20 .32 .32
18 cesd18 .47 .33 .69 -.40 .43 .66 .56 -.41 .57 .57
19 cesd19 .24 .17 .35 -.31 .32 .35 .32 -.24 .48 .37
20 cesd20 .35 .32 .52 -.30 .51 .55 .51 -.25 .38 .32
cesd11 cesd12 cesd13 cesd14 cesd15 cesd16 cesd17 cesd18 cesd19 cesd20
1 .28 -.41 .33 .26 .23 -.31 .27 .47 .24 .35
2 .25 -.35 .28 .31 .09 -.33 .29 .33 .17 .32
3 .34 -.64 .40 .45 .16 -.60 .49 .69 .35 .52
4 -.18 .49 -.23 -.32 -.09 .48 -.18 -.40 -.31 -.30
5 .21 -.38 .29 .26 .10 -.31 .30 .43 .32 .51
6 .37 -.65 .35 .49 .16 -.56 .36 .66 .35 .55
7 .30 -.53 .29 .39 .21 -.43 .37 .56 .32 .51
8 -.19 .56 -.21 -.34 -.07 .53 -.20 -.41 -.24 -.25
9 .21 -.43 .26 .52 .28 -.40 .32 .57 .48 .38
10 .30 -.44 .33 .38 .17 -.35 .32 .57 .37 .32
11 -.33 .19 .21 .16 -.24 .24 .31 .20 .35
12 -.33 -.39 -.46 -.17 .75 -.34 -.63 -.32 -.45
13 .19 -.39 .37 .15 -.32 .28 .40 .26 .33
14 .21 -.46 .37 .26 -.43 .36 .59 .44 .41
15 .16 -.17 .15 .26 -.15 .12 .19 .41 .16
16 -.24 .75 -.32 -.43 -.15 -.29 -.54 -.32 -.43
17 .24 -.34 .28 .36 .12 -.29 .57 .32 .30
18 .31 -.63 .40 .59 .19 -.54 .57 .40 .47
19 .20 -.32 .26 .44 .41 -.32 .32 .40 .36
20 .35 -.45 .33 .41 .16 -.43 .30 .47 .36
→再現性問題や統計の誤使用の問題からデータ収集中の検定の実施は行うべきではない(有意になるまでNを増やすN増しは絶対にダメ)。データ収集前に必要なサンプルサイズを決めておく必要がある。
→検定力が低い研究は,帰無仮説の棄却ができず結論づけができないので,研究を行う意義が薄くなる。
必要なサンプルサイズは,効果量,検定力,有意水準の3つを設定することで計算できる。
検定力と有意水準は,慣例に従い,検定力を0.8,有意水準を0.05に設定することが多い。効果量については,研究の対象やアウトカム(効果指標)に応じて設定していく必要がある(文脈によっては小さい効果量に意味があることもある)。
→効果量を決めることが一番むずかしい。効果量が決まれば,サンプルサイズも決まってくる。
Hislop et al. (2014)は,過去の無作為化比較試験で行われたサンプルサイズ設計を系統的に展望して,効果量の設定方法について整理した。
(1)治療者や患者が実際に差を感じるような意味づけのなされた差を採用する「重要な差(important difference)」(アンカーに基づく方法,分布に基づく方法,医療経済学に基づく方法,標準化効果量に基づく方法)
(2)過去の研究における群間差を参考にする「現実的な差(realistic difference)」(パイロット研究)
(3)重要な差と現実的な差の両方を含む決定方法として,意見聴取に基づく方法とエビデンスの展望に基づく方法がある
→絶対的な基準はないので,研究領域で納得感のある決め方ができれば良い。
2群の平均値差の必要サンプルサイズ
library(pwr)
pwr.t.test(n = NULL,
d = 0.3,
sig.level = 0.05,
power = 0.8,
type = "two.sample",
alternative = "two.sided")
Two-sample t test power calculation
n = 175.3847
d = 0.3
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
library(pwr)
pwr.t.test(n = NULL,
d = 0.6,
sig.level = 0.05,
power = 0.8,
type = "two.sample",
alternative = "two.sided")
Two-sample t test power calculation
n = 44.58577
d = 0.6
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
相関分析の必要サンプルサイズ
approximate correlation power calculation (arctangh transformation)
n = 193.0867
r = 0.2
sig.level = 0.05
power = 0.8
alternative = two.sided
approximate correlation power calculation (arctangh transformation)
n = 45.91614
r = 0.4
sig.level = 0.05
power = 0.8
alternative = two.sided
脱落を考慮して調査を計画する
必要サンプルサイズは,帰無仮説検定を行う上で必要なサンプルサイズなので,データ収集完了時にこれを下回らないようにする必要がある。
脱落を考慮して調査計画を立てれば安心して実施できる。そんなに難しい話ではなく,暫定的に脱落率を決めて(オンライン調査だと30%など),必要サンプルサイズに\(\frac{1}{1-脱落率}\) をかければ良い。
66名からデータを得るように計画すれば,3割が脱落しても45名以上からデータが得られる。
Woodworth et al.(2018)のデータについて,各種検定を実施してみよう!ご自身の研究についてサンプルサイズ設計をしてみよう!
回帰でgtsummaryを使うと便利だね。
mod1 <- glm(response ~ trt + age + grade, trial, family = binomial)
t1 <- tbl_regression(mod1, exponentiate = TRUE)
1.モデルの特定
→結果変数と説明変数を連結する関係式&反応変数の確率分布の2つを指定
2.モデルのパラメータの推定
3.モデルの妥当性のチェック
→モデルはどのくらいよくデータに当てはまり、データを要約しているか
4.推測
→モデルのパラメータに関する仮説検定と信頼区間の計算および結果の解釈
→独立性のないデータは、一般化線形混合モデルを用いる。
ヒストグラム、歪度、尖度、Kolmogorov-Smirnov(コロモゴロフ・スミノフ)検定などで確認
データの分布によっては対数変換、外れ値の吟味
→現実にあわないデータ処理によって、データが歪む可能性もあるので注意
平均値差の検定の場合は、Welch補正を使う。一般線形モデルでは、ロバスト回帰分析もある。
一般線形モデルの仮定が難しいデータの場合は、ノンパラメトリック分析へ
反復測定などにより階層性のあるデータは、一般(化)線形混合モデル(Generalized Linear Mixed Model: GLMM)を用いる。
一般(化)線形混合モデルでは、固定効果(fixed effects)とランダム効果(random effects)の両方が含まれる。
集団全体の傾きや切片(固定効果)に、切片と傾きの個人差(ランダム効果)を考慮するモデルになる。
ランダム傾きモデル:切片は全体と同じだが、傾きが個人によって異なる
ランダム切片&傾きモデル:切片も傾きも個人によって異なる。
ランダムな欠測には影響をうけない(利用可能なデータをすべて使える)
不等分散や非独立性に強い
少サンプルでもバイアスが少ない(正確さが増すので検出力も上がる)
査読で求められることが増えてきている・・・
はしらないようにしてある。
library(lme4)
library(phia)
library(lmerTest)
BtheBlong %>%
dplyr::group_by(treatment,time) %>%
dplyr::summarise(Mean=mean(bdi, na.rm=T), SD=sd(bdi, na.rm=T), Min=min(bdi, na.rm=T), Max=max(bdi, na.rm=T))
BtheBlong %>%
ggplot(aes(x=time, y=bdi, colour=treatment, group=id))+
geom_line(alpha=.5)+
scale_x_discrete(label=c("Post","1mFU","3mFU","6mFU"),name="Time")+
scale_y_continuous(name="BDI")
#ランダム切片モデルを使った一般化線形混合モデル解析を実施した(実際は複数のモデルで比較する必要がある)。
result <- lmer(bdi~bdi.pre + time*treatment+(1|id),data = BtheBlong, REML = FALSE,na.action = na.omit)
summary(result)
result95ci <- confint(result,method = "Wald")
result95ci
resultAdjM <- interactionMeans(result)
print(resultAdjM)
resultAdjM %>%
dplyr::rename(adjMean = `adjusted mean`, se = `SE of link`) %>%
ggplot(aes(x = time, y = adjMean, colour = treatment, group = treatment)) +
geom_line(alpha=.5)+theme_classic(base_size = 12,base_family = "") +
geom_errorbar(aes(ymin = adjMean - 1.96*se, ymax = adjMean + 1.96*se,width = 0.1))+
scale_x_discrete(label=c("Post","1mFU","3mFU","6mFU"),name="Time")+
scale_y_continuous(name="Adjusted Mean") +
labs(colour = "Treatment")構造方程式モデル
| 潜在変数 | ||
|---|---|---|
| 観測変数 | 量的変数 | 質的変数 |
| 量的変数 | 因子分析 | 潜在プロフィール分析 |
| 質的変数 | 項目反応理論 | 潜在クラス分析 |
媒介分析
確証的因子分析
奥村さんのコード
#----------------------------#
#á@ÉfÅ[É^ÇÃì«Ç›çûÇ›
#----------------------------#
##(1) égópÉâÉCÉuÉâÉäÇÃèÄîı
install.packages(c("psych", "paran"), contriburl=contrib.url("http://cran.md.tsukuba.ac.jp/")) #1âÒÉCÉìÉXÉgÅ[ÉãÇ∑ÇÍÇŒ2ìxñ⁄ÇÕïsóv
library(psych)
library(paran)
##(2) ì«Ç›çûÇ›
setwd("W:/workshop2011/archives/04/")
dat <- read.table("HolzingerGW.csv", sep=",", header=T)
head(dat)
tail(dat)
#----------------------------#
#áAàÍïœó âêÕÇ∆ìÒïœó âêÕ
#----------------------------#
summary(dat)
describe(dat)
cor(dat)
boxplot(dat)
graphics.off()
#----------------------------#
#áBàˆéqêîÇÃåàíË
#----------------------------#
##(1) ÉJÉCÉUÅ[äÓèÄÇ∆ÉXÉNÉäÅ[äÓèÄ
x <- eigen(cor(dat))
plot(x$values)
abline(h=1)
graphics.off()
##(2) ïΩçsï™êÕ
paran(dat, cfa=TRUE, all=TRUE, graph=TRUE, iterations=1000)
graphics.off()
##(3) ç≈è¨ïΩãœïŒëää÷
VSS(dat, n= 4, plot=FALSE)
##(4) BICÇ∆RMSEA
fit1 <- fa(r=dat, nfactors=1, fm="ml", rotate="promax")
fit2 <- fa(r=dat, nfactors=2, fm="ml", rotate="promax")
fit3 <- fa(r=dat, nfactors=3, fm="ml", rotate="promax")
rbind(fit1$RMSEA, fit2$RMSEA, fit3$RMSEA)
rbind(fit1$BIC, fit2$BIC, fit3$BIC)
##(5) ã§í ê´ÅCàˆéqç\ë¢ÅCàˆéqÉpÉ^ÉìÅCàˆéqä‘ëää÷
fit2
fit2$communality #ã§í ê´
fit2$Structure #àˆéqç\ë¢
print(fit2$loadings, cutoff=0) #àˆéqÉpÉ^Éì
fit2$Phi #àˆéqä‘ëää÷系統的展望のポイントについては,1日目の批判的吟味を参照ください。ここでは,解析部分のメタ分析に焦点を当てます。
個々の研究の効果量を計算
個々の研究の効果量を統合
統計学的異質性を評価
事前に計画したサブグループ解析
感度分析
非報告バイアス(出版バイアス)の検討
2群の差は,標準化平均値差(Hedgesのg)を用いる
\(Hedgesのg = \frac{介入群の平均-統制群の平均}{2群をプールした標準偏差}\)
\(2群をプールした標準偏差=\sqrt{\frac{(n_{介入}-1)\hat{\sigma}^{2}_{介入}+(n_{統制}-1)\hat{\sigma}^{2}_{統制}}{(n_{介入}+n_{統制}-2)}}\)
\(\hat{\sigma}^{2}\)は不偏分散。他にCohenのdやGlassの⊿もあるが,バイアスを修正した効果量のHedgesのgがよく使用される
\(標準化平均値差の分散=\frac{n_{介入}+n_{統制}}{n_{介入}n_{統制}}+\frac{g^2}{2(n_{介入}+n_{統制})}\)
\(標準化平均値差の標準誤差=\sqrt{標準化平均値差の分散}\)
分散と標準偏差は効果量とサンプルサイズで決まる
2変数間の関連の強さは,ピアソンの積率相関係数を用いる
\(r = \frac{共分散_{x,y}}{標準偏差_{x}標準偏差_{y}}\)
\(相関係数の分散 = \frac{(1-r^2)^2}{n-1}\),\(相関係数のSE = \frac{1-r^2}{\sqrt{n-1}}\)
統計学的特性の望ましさから(サンプルサイズが小さいとrのseにはバイアスが生じる),フィッシャーのZ変換を用いることが多い
\(z = 0.5ln(\frac{1+r}{1-r})\), \(zの分散 = \frac{1}{n-3}\), \(zのSE = \frac{1}{\sqrt{n-3}}\)
zからrへは\(\frac{exp(2z)-1}{exp(2z)+1}\)で変換できる(信頼区間も同様)。
リスク差=介入群のありの比率-統制群のありの比率
\(治療必要数=\frac{1}{リスク差}\)
\(リスク比=\frac{介入群のありの比率}{統制群のありの比率}\)
\(オッズ比=\frac{\frac{介入群のありの比率}{1-介入群のありの比率}}{\frac{統制群のありの比率}{1-統制群のありの比率}}\)
リスク差や治療必要数(Number needed to treat)は解釈しやすいが,計算上オッズ比が便利なので,オッズが使用される。さらに対数オッズ比(ln(オッズ比))を使う(1を起点に解釈→0を起点に解釈,分散も計算がしやすい)。
オッズ比の分散=介入群あり比率+介入群なし比率+統制群あり比率+統制群なし比率
1つの真の効果を仮定し,個々の研究の効果は真の効果と研究内分散(標本誤差)によるとするモデル
General variance based法による平均値差の統合の場合,重みは付けは,各研究の効果量の分散の逆数を使う(効果量の分散はサンプルサイズに依存)。
\(統合された効果量=\frac{\Sigma W_{i}Y_{i}}{\Sigma W_{i}}=\frac{(各研究の重みW_{i}*各研究の効果量Y_{i})の合計}{各研究の重みW_{i}の合計}\)
\(W_{i} = \frac{1}{V_{i}} = \frac{1}{効果量の分散V_{i}}\)
\(統合された効果量の標準誤差=\sqrt{\frac{1}{各研究の重みW_{i}の合計}}\)
\(95\%CIの下限=統合された効果量-1.96*統合された効果量の標準誤差\)
\(95\%CIの上限=統合された効果量+1.96*統合された効果量の標準誤差\)
真の効果を分布として考えて,研究ごとに真の効果が異なると仮定する。
個々の研究の効果は,真の効果の分布の平均値(μ)と研究間分散( \(\tau^2\))と研究内分散(標本誤差)によるとするモデル
1.固定効果モデルを用いて統合された効果量を算出し,Qを計算する(Tは個々の研究の効果量\(\bar{T}\)は統合された効果量)。Qは研究間のばらつきの程度。
\(Q = \Sigma W_{i}(T_{i}-\bar{T})^2\)
2.Qを用いて,研究間分散(\(\tau^2\))を計算する(kはメタ分析に含めた研究数)
\(\tau^2 = \frac{Q-(k-1)}{\Sigma W_{i} - \frac{\Sigma W^2_{i}}{\Sigma W_{i}}}\)
3.重み付けの計算に\(\tau^2\)\(V_{i}\)を用いる以外は,以降は固定効果モデルと同じ。
\(W_{i} = \frac{1}{V_{i}+\tau^2}\)
DerSimonian-Lairdの他に,Paule-Mandel, Restricted Maximum-Likelihood, Maximum-likelihood, Hunter-Schmidt, Sidik-Jonkman, Hedges, Empirical Bayesがある(Rのmeta パッケージで利用可能)。
DerSimonian-Lairdがよく使われるが,それはRevManで使われているなどの理由であって,Maximum-Likelihood, Sidik-Jonkman, Empirical Bayesの方が研究間変動の推定には良い性質がある(Harrer rt al., 2019)
ランダム効果か固定効果かは事前登録の段階で決める。
概念的異質性
論文を精査することで明らかとなる質的な研究間の違い
臨床的異質性:年齢,性別,人種,重症度,合併症など患者集団の違い(心理学だと対象者属性の異質性)
方法論的異質性:治療法,治療期間,治療者の経験など方法の違い(心理学だと測定法,実験操作などの違い)
統計学的異質性
→統計学的異質性も概念的異質性もある場合:(1)メタ分析しない,(2)サブグループ解析・メタ回帰分析の実施
→ 検定ではなく,Higgins et al.(2003)の\(I^2=\frac{Q-(k-1)}{Q}100\)のような記述的指標を用いて判断をする。 \(I^2\)は,効果量のバラツキの内,効果量の異質性が占める割合を%で表している。
\(I^2\)のおおまかな目安
40%以下 : 異質性は問題ではないかも
30~60% : 中程度の異質性があるかも
50~90% : 実際に異質性があるかも
75~100% : 無視できない異質性がある
→プロトコルの段階で計画を立てる
感度分析
→(1)適格基準の変更(患者の特徴を変えても変わらない?),(2)分析データの変更(調整済み効果量では?),(3)分析法の変更(固定効果とランダム効果モデルで大きな差はない?)
ネガティブな結果が出版されにくいと,効果を過剰評価するかもしれない(Funnel plotやEggerの検定などで出版バイアスを確認)。
出版バイアスがある場合は,Trim-fill法などで効果の過大評価を調整する。
Trim-fill法:偏った側からデータを除外して,左右が対称になるように調整する。その後,除外したデータを戻して,その反対側にデータを追加する
メタ分析をSEMの枠組みで実行することができる。
SEMを活用することで,多変量メタ分析ができる。例えば,臨床試験には複数のアウトカムが設定されているが,アウトカム間に関連があることもあり,それを考慮した上で効果を検討することができる。
メタ分析によって複数の研究で報告されている相関行列を統合し,SEMで確証的因子分析を実施したり,媒介分析を実施したりもでる。
詳しくは,Jak(2015),Cheung(2015),宇佐美(2018)を参照
RのmetaSEMパッケージが使える。
# パッケージの読みこみ
library(dmetar)
library(meta)
# データを作る(Mは平均,Sは標準偏差,Nは人数,eが介入・実験群,cが統制群)
Me <- c(35, 40, 35, 22, 38, 31, 29, 28)
Se <- c(7.5, 8.2, 5.9, 10.4, 8.8, 9.1, 6.4, 4.6)
Ne <- c(15, 32, 18, 22, 27, 42, 29, 17)
Mc <- c(20, 23, 30, 10, 22, 25, 19, 20)
Sc <- c(6.5, 4.3, 5.7, 7.8, 6.9, 10.4, 7.1, 6.5)
Nc <- c(16, 30, 19, 21, 29, 40, 32, 18)
Author<- c("Senshu (1923)","Waseda (1920)","Doshisha (1920)","Kangaku (1932)","Keio (1920)","Ritsumei (1922)","Hyokyo (1978)","Toyo (1928)")
Subgroup <- c("East","East","West","West","East","West","West","East")
# データフレームを作る
data_meta <- data.frame(Author, Me, Se, Ne, Mc, Sc, Nc, Subgroup)metaパッケージの連続変数用メタ分析関数のmetacontを使う。
comb.fixedが固定効果,comb.randomがランダム効果になる。
method.smdは,デフォルトが”Hedges”だが,“Cohen”や”Glass”も選べる。
method.tauは,デフォルトがDL(DerSimonian-Laird)だが他にも選択できる。今回は,最尤法(ML)を選択した。
Number of studies: k = 8
Number of observations: o = 407
SMD 95%-CI z p-value
Common effect model 1.3693 [1.1470; 1.5917] 12.07 < 0.0001
Random effects model 1.4897 [1.0553; 1.9240] 6.72 < 0.0001
Quantifying heterogeneity:
tau^2 = 0.2783 [0.0839; 1.6462]; tau = 0.5276 [0.2897; 1.2831]
I^2 = 77.7% [55.8%; 88.7%]; H = 2.12 [1.50; 2.97]
Test of heterogeneity:
Q d.f. p-value
31.33 7 < 0.0001
Details on meta-analytical method:
- Inverse variance method
- Maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Hedges' g (bias corrected standardised mean difference; using exact formulae)
SMD 95%-CI %W(common) %W(random)
Senshu (1923) 2.0867 [1.1906; 2.9828] 6.2 10.1
Waseda (1920) 2.5400 [1.8611; 3.2188] 10.7 12.3
Doshisha (1920) 0.8437 [0.1679; 1.5195] 10.8 12.4
Kangaku (1932) 1.2770 [0.6157; 1.9383] 11.3 12.5
Keio (1920) 2.0041 [1.3542; 2.6540] 11.7 12.6
Ritsumei (1922) 0.6093 [0.1658; 1.0527] 25.1 14.9
Hyokyo (1978) 1.4568 [0.8878; 2.0257] 15.3 13.5
Toyo (1928) 1.3813 [0.6352; 2.1273] 8.9 11.6
Number of studies: k = 8
Number of observations: o = 407
SMD 95%-CI z p-value
Common effect model 1.3693 [1.1470; 1.5917] 12.07 < 0.0001
Random effects model 1.4897 [1.0553; 1.9240] 6.72 < 0.0001
Quantifying heterogeneity:
tau^2 = 0.2783 [0.0839; 1.6462]; tau = 0.5276 [0.2897; 1.2831]
I^2 = 77.7% [55.8%; 88.7%]; H = 2.12 [1.50; 2.97]
Test of heterogeneity:
Q d.f. p-value
31.33 7 < 0.0001
Details on meta-analytical method:
- Inverse variance method
- Maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Hedges' g (bias corrected standardised mean difference; using exact formulae)
異質性に関連した外れ値の検討,Influence Analyses,GOSH Plot Analysisは省略しました。詳細は,Doing Meta-Analysis in R : A Hands-on Guideを確認ください
dmetarパッケージのeggers.testを使う。
Number of studies: k = 10 (with 2 added studies)
Number of observations: o = 500
SMD 95%-CI z p-value
Random effects model 1.2193 [0.7183; 1.7204] 4.77 < 0.0001
Quantifying heterogeneity:
tau^2 = 0.5286 [0.2202; 2.3463]; tau = 0.7271 [0.4692; 1.5318]
I^2 = 83.5% [71.1%; 90.5%]; H = 2.46 [1.86; 3.25]
Test of heterogeneity:
Q d.f. p-value
54.39 9 < 0.0001
Details on meta-analytical method:
- Inverse variance method
- Maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Trim-and-fill method to adjust for funnel plot asymmetry (L-estimator)
今回触れていないこと
サンプルサイズを考慮して真の効果量を検討するP-curve分析(Simonsohn et al., 2014)
Risk of Biasのプロット方法
サブグループ解析
2値変数,相関係数のメタ分析
ベイジアンメタ分析
ネットワークメタ分析
メルチレベルメタ分析
Doing Meta-Analysis in R : A Hands-on Guideでより詳しく学ぶことができます。
関心のあるテーマでオープンデータ・模擬データを探して,Quartoを使って解析レポートを作成しましょう。